

    ##############################################################################################
    ### To produce Figure 4 (Senate) and Figure A1: requires "count1_sen.csv" and "count_sen.bug" 
    ##############################################################################################

    # needs these packages
    library(R2WinBUGS)
    library(coda)    
 
    setwd("C:\\temp1") # Make necessary change

    T1 <- 102
    T2 <- 114
    
    Micro_Data <-  read.csv("count1_sen.csv", header=T)
    
    GOP_prez <- c(1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0) # from 102th to 114th Congresses

   
    ### Setup data

    # Dependent Variable    
    Y   <- Micro_Data$count  
     
    # Member-level IVs
    GOP <- ifelse(Micro_Data$party==200, 1, 0)
    seniority <- Micro_Data$Seniority_Year
    seniority.std <- as.vector((seniority - mean(seniority))/sd(seniority))
    fb <- Micro_Data$foreign_born2005/Micro_Data$pop2005
    fb.std <- as.vector((fb - mean(fb))/sd(fb))
    veteran <- Micro_Data$veteran
    Southern.State <- c(40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 51, 54) # VA, AL, AR, FL, GA, LA, MS, NC, SC, TX, plus KY and TN
    South <- ifelse(is.element(Micro_Data$ICPSR_state, Southern.State), 1, 0)
   
    # For the Poisson process
    X1 <- GOP
    X2 <- Micro_Data$FR_cmt
    X3 <- Micro_Data$AS_cmt
    X4 <- seniority.std
    X5 <- fb.std
    X6 <- veteran
    X7 <- South

    beta.label <-  c("SFRC", "SASC", "Seniority", "Foreign-Born(%)", "Veteran", "South")
    n.beta  <- 6

    congress <- Micro_Data$congress - (min(Micro_Data$congress) - 1)
    alpha.label <- paste(T1:T2, "", sep="")
    zeta.label <- paste(T1:T2, "", sep="")
    n.zeta <- length(zeta.label)

    # For the Binary process
    K1 <- Micro_Data$FR_cmt

    a.label <- c("Intercept", "SFRC")
    n.a     <- 2

    # Congress-level IVs
    Z1 <- GOP_prez 

    gamma.label <- c("Intercept", "GOP Presidency")
    n.gamma <- 2
    zeta.gamma.label <- c("Intercept", "GOP Presidency")
    n.zeta.gamma <- 2

    N <- length(Y)
    J <- length(Z1)
   
    data <- list("N", "J", "n.beta","n.gamma","n.zeta.gamma", "n.a", "Y", 
                 "X1", "X2", "X3", "X4", "X5", "X6", "X7", "congress", "K1",  
                "Z1")

    # setup initial values
    init1 <- list(beta=rnorm(n.beta), 
                  a=rnorm(n.a), zeta=rnorm(J), zeta.gamma=rnorm(n.zeta.gamma),
                  gamma=rnorm(n.gamma), sigma.zeta=runif(1),
                  alpha=rnorm(J), sigma.alpha=runif(1), u=sample(c(0,1), N, replace=T)) 
    init2 <- list(beta=rnorm(n.beta), 
                  a=rnorm(n.a), zeta=rnorm(J), zeta.gamma=rnorm(n.zeta.gamma),
                  gamma=rnorm(n.gamma), sigma.zeta=runif(1),
                  alpha=rnorm(J), sigma.alpha=runif(1), u=sample(c(0,1), N, replace=T))    
    init3 <- list(beta=rnorm(n.beta), 
                  a=rnorm(n.a), zeta=rnorm(J), zeta.gamma=rnorm(n.zeta.gamma),
                  gamma=rnorm(n.gamma), sigma.zeta=runif(1),
                  alpha=rnorm(J), sigma.alpha=runif(1), u=sample(c(0,1), N, replace=T))

    inits <- list(init1, init2, init3)
    
    # setup parameters
    parameters = c("beta", "a", "zeta", "zeta.gamma", "gamma", "alpha", "sigma.zeta", "sigma.alpha")
    
    # Run bugs.  
    Model.fit = bugs(data, inits, parameters, "count_sen.bug", working.directory=getwd(), 
                        n.chains=3, n.thin=5, n.burnin=2000, n.iter=7000)

    Rhats <- Model.fit$summary[,8]
    write.table(Rhats, "Rhats.txt")
    range(Rhats)
    

    ### 4. Analyze the draws from the Posterior distribution ############################################

    chain1 <- read.coda(output.file="coda1.txt", index.file="codaIndex.txt", quiet=T)
    chain2 <- read.coda(output.file="coda2.txt", index.file="codaIndex.txt", quiet=T)
    chain3 <- read.coda(output.file="coda3.txt", index.file="codaIndex.txt", quiet=T)

    MCMC <- rbind(chain1, chain2, chain3)
    iter <- nrow(MCMC)

    a.finder <- seq(from=1, to=n.a, by=1)
    alpha.finder  <- seq(from=(n.a+1), to=(n.a+J), by=1)
    beta.finder <- seq(from=(n.a+J+1), to=(n.a+J+n.beta), by=1)
    gamma.finder  <- seq(from=(n.a+J+n.beta+2), to=(n.a+J+n.beta+1+n.gamma), by=1)
    zeta.finder  <- seq(from=(n.a+J+n.beta+1+n.gamma+3), to=(n.a+J+n.beta+1+n.gamma+2+J), by=1)
    zeta.gamma.finder <- seq(from=(n.a+J+n.beta+1+n.gamma+2+J+1), to=(n.a+J+n.beta+1+n.gamma+2+J+ n.zeta.gamma), by=1)

    a.mc  <- MCMC[,a.finder]
    alpha.mc  <- MCMC[,alpha.finder]
    beta.mc   <- MCMC[,beta.finder]
    gamma.mc  <- MCMC[,gamma.finder]
    zeta.mc   <- MCMC[,zeta.finder]
    zeta.gamma.mc  <- MCMC[,zeta.gamma.finder]

    ### Credible Intervals
    a.CI      <- round(apply(a.mc,      2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    beta.CI   <- round(apply(beta.mc,   2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    gamma.CI  <- round(apply(gamma.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    alpha.CI  <- round(apply(alpha.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
    zeta.CI   <- round(apply(zeta.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995, 0.1, 0.9)), 3)
    zeta.gamma.CI  <- round(apply(zeta.gamma.mc,  2, FUN=quantile, probs=c(0.025, 0.975, 0.5, 0.05, 0.95, 0.005, 0.995)), 3)
 

    pdf("figure4_senate.pdf", height=8, width=5)
    par(mfrow=c(1,1), mar = c(0,0,0,0), oma = c(3,4,0,0), mgp=c(1.5, 0, 0), xpd=NA, tck=-0.007)  #c(bottom, left, top, right)
    plot(1:n.zeta, 1:n.zeta, type='n', col='black', xlim=c(-1, 1.5), ylim=c(0,n.zeta+1),
           xlab="", ylab="", yaxt="n") # , yaxt="n"
    text(    zeta.CI[3,], 1:n.zeta-0.25, labels=zeta.label,  cex=1) # use srt for text rotation 90=degree
    points(  zeta.CI[3,1:n.zeta], 1:n.zeta, type='p', pch=19, col='blue') # Plot the Medians
    segments(zeta.CI[1,1:n.zeta], 1:n.zeta, 
             zeta.CI[2,1:n.zeta], 1:n.zeta, col='red', lwd=2) # Plot CI
    segments(zeta.CI[4,1:n.zeta], 1:n.zeta, 
             zeta.CI[5,1:n.zeta], 1:n.zeta, col='black', lwd=3) # Plot CI
    legend('top', legend="Varying-Slopes of GOP")
    #points(  rep(0, n.zeta), 1:n.zeta, type='l', lty=2) # use this instead of abline 
    abline(v=c(0), lty=2)    
    segments(0, 0.5, 0, length(zeta.label)+0.5, lty=2)
    segments(-1, 1.5, 1.5, 1.5, lty=2) # Bush I
    text(-0.7, 0.8, labels="Bush I", srt=90)
    segments(-1, 5.5, 1.5, 5.5, lty=2) # Clinton
    text(-0.7, 3.5, labels="Clinton", srt=90)
    segments(-1, 9.5, 1.5, 9.5, lty=2) # Bush II
    text(-0.7, 7.5, labels="Bush II", srt=90)
    text(-0.7, 11.5, labels="Obama", srt=90)
    dev.off() 


    plan <- rbind(c(1, 5), c(2, 5), c(3, 5), c(4, 5))
    print(plan)

    pdf("figure_A1.pdf", height=7, width=5)
    layout(plan)
    par(mar = c(2,1,2.5,1), oma = c(3,4,0,0), mgp=c(1.5, 0, 0), xpd=NA, tck=-0.007)  #c(bottom, left, top, right)

    plot(1:n.a, 1:n.a, type='n', col='black', xlim=c(-1.5, 1.5), ylim=c(0.5,n.a+0.5),
            xlab="", ylab="", yaxt="n")
    text(    a.CI[3,1:n.a], 1:n.a-0.25, labels=a.label,  cex=1) # use srt for text rotation 90=degree
    points(  a.CI[3,1:n.a], 1:n.a, type='p', pch=19, col='blue') # Plot the Medians
    segments(a.CI[1,1:n.a], 1:n.a,  
             a.CI[2,1:n.a], 1:n.a, col='red', lwd=2) # Plot 95% CI 
    segments(a.CI[4,1:n.a], 1:n.a,  
             a.CI[5,1:n.a], 1:n.a, col='black', lwd=3) # Plot 90% CI
    title(main=expression(paste("(1) Binary Process (", omega, ")", sep="")), font.main=1)
    #title(main="(1) Binary Process", font.main=1)
    segments(0, 0.5, 0, length(a.label)+0.5, lty=2)
    
    plot(1:n.beta, 1:n.beta, type='n', col='black', xlim=c(-1.5, 1.5), ylim=c(0.5,n.beta+0.5), 
            xlab="", ylab="", yaxt="n")
    text(    beta.CI[3,1:n.beta], 1:n.beta-0.4,  labels=beta.label,  cex=1) # use srt for text rotation 90=degree
    points(  beta.CI[3,1:n.beta], 1:n.beta, type='p', pch=19, col='blue') # Plot the Medians
    segments(beta.CI[1,1:n.beta], 1:n.beta,  
             beta.CI[2,1:n.beta], 1:n.beta, col='red', lwd=2) # Plot 95% CI 
    segments(beta.CI[4,1:n.beta], 1:n.beta,  
             beta.CI[5,1:n.beta], 1:n.beta, col='black', lwd=3) # Plot 90% CI
    title(main=expression(paste("(2) Poisson Process (", beta[2]," - ", beta[7], ")", sep="")), font.main=1)
    #title(main="(2) Poisson Process", font.main=1)
    segments(0, 0.5, 0, length(beta.label)+0.5, lty=2)

    plot(1:n.gamma, 1:n.gamma, type='n', col='black', xlim=c(-1.5, 1.5), ylim=c(0.5,n.gamma+0.5),
           xlab="", ylab="", yaxt="n") # , yaxt="n"
    text(    gamma.CI[3,], 1:n.gamma-0.25, labels=gamma.label,  cex=1) # use srt for text rotation 90=degree
    points(  gamma.CI[3,1:n.gamma], 1:n.gamma, type='p', pch=19, col='blue') # Plot the Medians
    segments(gamma.CI[1,1:n.gamma], 1:n.gamma, 
             gamma.CI[2,1:n.gamma], 1:n.gamma, col='red', lwd=2) # Plot CI
    segments(gamma.CI[4,1:n.gamma], 1:n.gamma, 
             gamma.CI[5,1:n.gamma], 1:n.gamma, col='black', lwd=3) # Plot CI
    title(main=expression(paste("(3) Predictors for Varying Intercepts (", gamma, ")", sep="")), font.main=1)
    #title(main="(3) Congress-Level Factors \n for Varying-Intercepts", font.main=1)
    segments(0, 0.5, 0, length(gamma.label)+0.5, lty=2)


    plot(1:n.zeta.gamma, 1:n.zeta.gamma, type='n', col='black', xlim=c(-1.5, 1.5), ylim=c(0.5,n.zeta.gamma+0.5),
           xlab="", ylab="", yaxt="n") # , yaxt="n"
    text(    zeta.gamma.CI[3,], 1:n.zeta.gamma-0.25, labels=zeta.gamma.label,  cex=1) # use srt for text rotation 90=degree
    points(  zeta.gamma.CI[3,1:n.zeta.gamma], 1:n.zeta.gamma, type='p', pch=19, col='blue') # Plot the Medians
    segments(zeta.gamma.CI[1,1:n.zeta.gamma], 1:n.zeta.gamma, 
             zeta.gamma.CI[2,1:n.zeta.gamma], 1:n.zeta.gamma, col='red', lwd=2) # Plot CI
    segments(zeta.gamma.CI[4,1:n.zeta.gamma], 1:n.zeta.gamma, 
             zeta.gamma.CI[5,1:n.zeta.gamma], 1:n.zeta.gamma, col='black', lwd=3) # Plot CI
    title(main=expression(paste("(4) Predictors for Varying Slopes (", delta, ")", sep="")), font.main=1)
    #title(mai="(4) Congress-Level Factors \n for Varying-Slope of GOP", font.main=1)
    segments(0, 0.5, 0, length(zeta.gamma.label)+0.5, lty=2)

    plot(1:n.zeta, 1:n.zeta, type='n', col='black', xlim=c(-1, 1.5), ylim=c(0.5,n.zeta+0.5),
           xlab="", ylab="", yaxt="n") # , yaxt="n"
    text(    zeta.CI[3,], 1:n.zeta-0.25, labels=zeta.label,  cex=1) # use srt for text rotation 90=degree
    points(  zeta.CI[3,1:n.zeta], 1:n.zeta, type='p', pch=19, col='blue') # Plot the Medians
    segments(zeta.CI[1,1:n.zeta], 1:n.zeta, 
             zeta.CI[2,1:n.zeta], 1:n.zeta, col='red', lwd=2) # Plot CI
    segments(zeta.CI[4,1:n.zeta], 1:n.zeta, 
             zeta.CI[5,1:n.zeta], 1:n.zeta, col='black', lwd=3) # Plot CI
    title(main=expression(paste("(5) Varying Slopes of GOP (", beta[1*j], ")", sep="")), font.main=1)
    #title(main="(5) Varying-Slopes of Republicans", font.main=1)
    segments(0, 0.5, 0, length(zeta.label)+0.5, lty=2)
    segments(-1, 1.5, 1.5, 1.5, lty=3) # Bush I
    text(-0.7, 1, labels="Bush I", srt=90)
    segments(-1, 5.5, 1.5, 5.5, lty=3) # Clinton
    text(-0.7, 3.5, labels="Clinton", srt=90)
    segments(-1, 9.5, 1.5, 9.5, lty=3) # Bush II
    text(-0.7, 7.5, labels="Bush II", srt=90)
    text(-0.7, 11.5, labels="Obama", srt=90)

    dev.off()
    
